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Abstract 

Bardet-Biedl syndrome (BBS) is a model ciliopathy characterized by a wide 
range of clinical variability. The heterogeneity of this condition is reflected in 
the number of underlying gene defects and the epistatic interactions between 
the proteins encoded. BBS is generally inherited in an autosomal recessive trait. 
However, in some families, mutations across different loci interact to modulate 
the expressivity of the phenotype. In order to investigate the magnitude of epis- 
tasis in one BBS family with remarkable intrafamilial phenotypic variability, we 
designed an exome sequencing-based approach using SOLID 5500x1 platform. 
This strategy allowed the reliable detection of the primary causal mutations in 
our family consisting of two novel compound heterozygous mutations in 
McKusick-Kaufman syndrome (MKKS) gene (p.D90G and p.V396F). Addition- 
ally, exome sequencing enabled the detection of one novel heterozygous NPHP4 
variant which is predicted to activate a cryptic acceptor splice site and is only 
present in the most severely affected patient. Here, we provide an exome 
sequencing analysis of a BBS family and show the potential utility of this tool, 
in combination with network analysis, to detect disease-causing mutations and 
second-site modifiers. Our data demonstrate how next-generation sequencing 
(NGS) can facilitate the dissection of epistatic phenomena, and shed light on 
the genetic basis of phenotypic variability. 
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Introduction 

Ciliopathies are an emerging group of clinical disorders 
characterized by large genetic heterogeneity and clinical 
overlap. Due to the ubiquitous nature of the primary cil- 
ium, ciliopathies can affect many organ systems (Fliegauf 
et al. 2007). Bardet-Biedl syndrome (BBS; MIM 209900) 
is considered a model of ciliopathy (Badano et al. 2006a), 
and its prevalence varies from 1:160,000 in Europe (Klein 
and Ammann 1969; Beales et al. 1997) to 1:13,500 and 
1:17,500 in Kuwait and Newfoundland, respectively (Farag 
and Teebi 1989; Green et al. 1989). BBS is a pleiotropic 
disorder that has variable expressivity and a wide range of 
clinical variability observed both within and between fam- 
ilies (Beales et al. 1999). Primary features include retinitis 
pigmentosa (RP), Polydactyly, obesity, genital defects, 
renal anomalies, and learning disabilities. Secondary fea- 
tures include speech delay, developmental delay, diabetes 
mellitus, dental anomalies, congenital heart disease, 
brachydactyly/syndactyly, ataxia/poor coordination, and 
anosmia/hiposmia (Forsythe and Beales 2013). 

The clinical diagnosis of BBS requires the presence of 
at least four primary or three primary and two secondary 
features (Beales et al. 1999; Rooryck and Lacombe 2008; 
Putoux et al. 2012). The spectrum of clinical features of 
BBS shares common characteristics with other ciliopathies 
such as Joubert syndrome (JBTS; MIM 213300), Leber 
congenital amaurosis (LCA; MIM 204000), McKusick- 
Kaufman syndrome (MKKS; MIM 236700), Meckel- 
Gruber syndrome (MKS; MIM 249000), nephronophthisis 
(NPH; MIM 256100), orofaciodigital syndrome type 1 
(OFD1; MIM 311200), and Senior-Loken syndrome (SLS; 
MIM 266900). Although these are considered as distinct 
clinical entities, in many cases it is very difficult to assign 
a specific clinical diagnosis due to phenotypic overlap. 
Accordingly, to the phenotypic overlap, genetic overlap 
also exists. Thus, mutations in MKKS (MIM 604896) also 
cause MKKS - characterized by genitourinary malforma- 
tions (hydrometrocolpos), Polydactyly, and more rarely 
gastrointestinal abnormalities (McKusick et al. 1964; 
Stone et al. 2000). 

To date, mutations in 18 genes have been associated 
with BBS (Retinal Information Network website. Avail- 
able: https://sph.uth.edu/retnet/) accounting for -70% of 
affected individuals (Muller et al. 2010). BBS is typically 
inherited in an autosomal recessive manner. However, the 
involvement of mutations in more than one locus has 
been substantially growing. Digenic triallelic inheritance 
has been reported in some patients with three mutations 
across two BBS loci that interact to cause disease (Katsan- 
is et al. 2001; Badano et al. 2003; Beales et al. 2003). In 
addition, the presence of second-site modifiers that may 
modulate the expression of the clinical phenotype has 



been proposed to explain, in part, the significant inter- 
and intrafamilial clinical variability in BBS (Katsanis 
2004; Badano et al. 2006b; Khanna et al. 2009). 

In the last few years, the introduction of the next-gener- 
ation sequencing (NGS) has revolutionized clinical genet- 
ics, making whole-exome sequencing (WES) a rapid way 
to elucidate the genetic basis of clinically and genetically 
heterogeneous Mendelian disorders (Bamshad et al. 2011; 
Ionita-Laza et al. 2011; Rabbani et al. 2012). At present, it 
is expected that NGS, in combination with network analy- 
sis (Minguez et al. 2009) and other advanced bioinformat- 
ics tools that allows prioritizing candidate genes because 
of their functional relationships (Ideker and Sharan 2008), 
will play an increasingly important role in the diagnosis of 
complex and oligogenic disorders. Moreover, WES pro- 
vides a unique possibility for investigating the presence of 
additional mutations that may modify the expressivity of 
BBS phenotype. Other authors have been interested in 
determining the contribution of epistasis in BBS (Katsanis 
et al. 2002; Hichri et al. 2005; Smaoui et al. 2006; Abu-Sa- 
fieh et al. 2012), but their approaches have been based on 
Sanger sequencing of known BBS loci (Katsanis et al. 
2002; Hichri et al. 2005; Smaoui et al. 2006; Abu-Safieh 
et al. 2012) or on targeted exon capture sequencing of a 
number of genes (Redin et al. 2012). 

In this study, we conducted a WES approach in a Span- 
ish family with significant intrafamilial variability in the 
BBS phenotype, from full blown to relatively mild forms 
of BBS. Using this technology, compound heterozygous 
mutations have been found in MKKS and the contribution 
of additional changes has been dissected. 

Material and Methods 

Subjects and clinical assessment 

Our study involved one Spanish family (RP42) com- 
prising three affected and four unaffected individuals 
(Fig. 1), all derived from the Ophthalmology Depart- 
ment to our Genetic Department. The study was carried 
out in accordance with the tenets of the Declaration of 
Helsinki and the ethical guidelines of our institutions. A 
group of 200 matching control individuals was also 
recruited. Written informed consent was obtained from 
all participants. 

Clinical diagnosis of retinal degeneration was based on 
visual acuity, fundus photography, computerized testing 
of central and peripheral visual fields, and electroretinog- 
raphy (ERG) findings. Typical ocular features include ini- 
tial night blindness, restriction of visual field, bone 
spicule pigmentation, attenuation of retinal vessels, waxy 
disk pallor, and abnormal ERG findings in a rod-cone 
pattern when recordable. 
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RP 42 



1:1 
[M2];[=] 
[M3];[=] 



l:2 
[M1];[=] 



11:1 
[M2];[=] 



ll:2 
[M1];[=] 
[M2];[=] 



□ 

ll:3 
[M1];[=] 



[M1];| 
[M2];| 
[M3];| 



ll:5 
[M1];[=] 
[M2];[=] 



M1: MKKS: c.269A>G p.D90G 
M2: MKKS: c.1 186G>T p.V396F 
M3: NPHP4: c.992+71G>T (splice site) 

Figure 1. Family cosegregation analysis. RP42 family tree showing 
the segregation of the sequence variants identified during the 
molecular analysis of MKKS and NPHP4. [M];[M]: homozygous; [M]; 
[=]: heterozygous. 

All subjects underwent a peripheral blood extraction 
for genomic DNA isolation from leukocytes using stan- 
dard protocols. DNA samples from individuals 1:1, 1:2, 
11:3, and 11:4 were processed for NGS. 

After the identification of the MKKS gene defects, 
affected individuals underwent clinical reevaluation 
focused on the identification of extraocular features asso- 
ciated with BBS. 

Description of DNA library preparation 
and sequencing 

Library preparation and exome capture were performed 
according to a protocol based on the Baylor College of 
Medicine protocol version 2.1 with several modifications. 
Briefly, 5 fig of input genomic DNA was sheared, end 
repaired, and ligated with specific adaptors. A fragment 
size distribution ranging from 160 bp to 180 bp after 
shearing and 200-250 bp after adaptor ligation was 
verified by Bioanalyzer (Agilent Technologies, Santa Clara, 
CA). The library was amplified by precapture linker- med- 
iated polymerase chain reaction (LM-PCR) using Fast- 
Start High Fidelity PCR System (Roche, Indianapolis, 
IN). After purification, 2 /ig of LM-PCR product was 
hybridized to NimbleGen SeqCap EZ Exome libraries V3. 
After washing, amplification was performed by postcapture 
LM-PCR using FastStart High Fidelity PCR System 
(Roche). Capture enrichment was measured by qPCR 
according to NimbleGen protocol. The successfully cap- 
tured DNA was measured by Quant-iT™ PicoGreen® 
dsDNA reagent (Invitrogen, Carlsbad, CA) and subjected 
to standard sample preparation procedures for sequenc- 
ing with SOLiD 5500x1 platform as recommended by the 



manufacturer. Shortly, emulsion PCR was performed on 
E80 scale (about 1 billion template beads) using a con- 
centration of 0.616 pm of enriched captured DNA. After 
breaking and enrichment, about 276 million enriched 
template beads were sequenced per lane on a six-lane 
SOLiD 5500x1 slide. 

Analysis of data from deep sequencing 

SOLiD 5500x1 reads were aligned against the human gen- 
ome reference (hgl9) using the program BFAST (Blat-like 
Fast Accurate Search Tool) allowing reads to map only to a 
unique position in the reference genome. Improperly 
mapped reads were filtered out with the SAMtools package, 
which was also used to sort SAM files and to generate and 
index BAM files. Variant calling was performed with the 
software GATK (Genome Analysis Toolkit) taking into 
account variants from NCBI database of Single Nucleotide 
Polymorphisms (dbSNP) for recalibration and realignment. 
Secondary analysis was performed by a custom shell script 
that queries VARIANT database (Medina et al. 2012), 
which includes SIFT (Kumar et al. 2009) and Polyphen-2 
(Adzhubei et al. 2010) scores, to annotate all single-nucleo- 
tide variations (SNVs) and small insertions and deletions 
(INDELS). Variants with an allele frequency higher than 
5% in the 1000 Genomes Project database were discarded. 
Only, exonic variants which produce a synonymous change 
in the open reading frame were discarded, whereas other 
type of variants in the exomes and variants in splice sites 
were kept for further analysis. Then, variants found in 
affected individuals were compared with variants present in 
not affected relatives. A last step was performed to compare 
the remaining variants with variants obtained from a group 
of healthy controls from the same local population as the 
family in study obtained from The Medical Genome Project 
(www.medicalgenomeproject.com). Finally, genes with 
variants in both alleles present in affected but not in healthy 
individuals of the family nor in the control local population 
and with a variant only in one allele in the carrier were 
ranked based on the analysis of the interactome, using Net- 
workMiner(Garcia-Alonso et al. 2012), to generate a list of 
candidate genes. 

Verification and assessment of the 
pathogenicity of variants 

Each predicted disease-causing variant was confirmed by 
Sanger sequencing, and cosegregation analysis was per- 
formed in the rest of the family members DNA samples. 

As mentioned above, we used Polyphen-2 and SIFT 
scores to evaluate the potential impact of novel missense 
substitutions on the function of the encoded protein. 
Evolutionary conservation across species was assessed 
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through the alignment of orthologous MKKS protein 
sequences (Pan troglodytes, Mus musculus, Canis lupus fa- 
miliaris, Callus gallus, Xenopus tropicalis, and Danio rerio) 
with the human MKKS protein sequence, using Clustal 
Omega Tool (Sievers et al. 2011). Furthermore, splice site 
tool Prediction by Neural Network (Reese 1997; Reese 
et al. 1997), exonic splicing enhancer prediction programs 
ESE Finder (Cartegni et al. 2003; Smith et al. 2006), and 
NetGene2 server (Hebsgaard et al. 1996) were applied to 
estimate the pathogenic nature of intronic sequence vari- 
ants that could affect the splicing process. The correct 
nomenclature for mutation was checked applying Muta- 
lyzer (Wildeman et al. 2008) using the corresponding 
Genbank reference sequences (MKKS; NG_009109.1 and 
NPHP4; NG_0 11724.2). Novel variants included in this 
article were submitted to the respective Locus Specific 
Database (LSDB) (http://grenada.lumc.nl/LOVD2/eye/home. 
php?select_db=MKKS). 

Network analysis 

Network enrichment analysis has been performed using 
the program SNOW (Minguez et al. 2009) included in 
the Babelomics (Medina et al. 2010) package (http:// 
www.babelomics.org). Genes are mapped onto the inter- 
actome (obtained from the STRING [Franceschini et al. 
2013] database), and the subnetwork connecting them is 
obtained. Several relevant parameters are calculated for 
this subnetwork, such as the connectivity or the number 
of components. An empirical distribution of the random 
expectation of these parameters is obtained by repeatedly 
sampling random sets of the same number of genes from 
the complete genome and calculating the average connec- 
tivity of their corresponding subnetworks. Thus, real val- 
ues of the parameters obtained for the genes analyzed can 
be contrasted with respect to their random expectations 
(Minguez et al. 2009; Minguez and Dopazo 2010). 

Results 

Clinical assessment 

The clinical findings of the affected individuals (11:2, 11:4, 
and 11:5) are reported in Table 1. Ocular manifestations 
of the disease in the three siblings were fairly typical of 
an early onset and severe form of RP. Night blindness 
was reported from the first decade. Thereafter, the disease 
rapidly progressed, and by the age of 20, 11:2, 11:4, and 
11:5 were severely visually disabled. The fundus examina- 
tion showed the typical signs of RP, including pale waxy 
disks, attenuation of the retinal vessels, and bone spicule 
pigmentation in the midperiphery. Photopic and scotopic 
ERGs were extinguished (11:2) or diminished (11:4) on 



Table 1. Summary of phenotype features documented in the 
affected members of the family RP 42. 



Feature 


ll:2 


Il:4 


!:5 


First symptom 


Night 


Night 


Night 




blindness 


blindness 


blindness 


Onset age 


3 


10 


8 


Primary features 








Retinitis pigmentosa 






+ 


Postaxial Polydactyly 

i J J J 


Both feet 


Both feet 


One foot 


\A/oinht nam snnmslu 
vvtriyriL uaii l aiiuilldiy 


Mnrmal 

I MUI I I Id I 




Di/onA/Din ht 
\J V tri Wtriu I I L 




weight 






Genital defects, 








hydrometrocolpos 








Learning disabilities 






+ 


Renal anomalies 








Kidney cysts 




- 




Kidney transplant 




- 




Secondary features 








Poor coordination 









Developmental delay 
Speech delay 
Brachydactyly/Syndactyly 
Dental anomalies 

Teeth crowding 

Congenital heart disease 



examination at ages 21 and 16, respectively. Patient 11:4 
showed extraocular features commonly associated with 
BBS such as postaxial Polydactyly, overweight, polycystic 
kidney, learning disabilities, and mild psychomotor delay. 
Patients 11:2 and 11:5 exhibited a less severe phenotype 
consisting of RP, postaxial Polydactyly, and mild learning 
disabilities. Subject 11:4 underwent renal transplantation 
due to polycystic kidney disease, whereas the other two 
did not manifest any renal abnormalities. Therefore, 
patient 11:4 fulfilled the criteria for the diagnosis of classic 
BBS, whereas 11:2 and 11:5 were diagnosed of mild BBS. 

Identification of causal mutations 

Upon exome sequencing, 43,570 sequence variants were 
identified in patient 11:4 (Table 2). The two most proba- 
bly pathogenic mutations located on exons 3 and 5 of 



Table 2. Variants identified 


by exome 


sequencing in the 


RP 42 


family. 










1:1 l:2 II : 3 ll:4 


Total SNVs 


45,068 


42,243 


43,269 


43,570 


Nonsynonymous SNVs 


6114 


5855 


5892 


5942 


Filtered dbSNP 


5979 


5735 


5756 


5821 


Filtered dbSNP and 1000 g 


5913 


5692 


5710 


5760 


Filtered dbSNP and 1000 g 


357 


373 


355 


367 


and predicted deleterious 











SNVs, single-nucleotide variations. 
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Chr20 



B 



Reverse strand 



MKKS 



29.04 Kb 



4 5 6 

H — ■ 



C Reference 



Sample 11:4 




Figure 2. Detection of novel mutations in the MKKS gene. (A) Chromosome overview of the chromosome 20, MKKS is mapped on region 
20p12.2 (black bar). (B) MKKS spans approximately 29 Kb and is composed of six exons. Filled boxes reflect coding exons (3-6) and unfilled 
boxes reflect UTR. (C) Electropherogram depiction of the index patient (ll:4) confirming the heterozygous mutations in exons 3 and 5 of MKKS 
gene. IUPAC SNP codes used to designate heterozygous substitutions ("k" from Keto for GfT and "r" from puRine for A/G). MKKS Genbank 
accession number: NG_009109.1. 



MKKS gene (Fig. 2) have never been reported in public 
variant databases such as dbSNP, Exome Variant Server 
(http://evs.gs.washington.edu/EVS/), or 1000 genomes 
(http://www.1000 genomes.org/). Sanger sequencing of 
these exons confirmed the novel missense variants 
(c.269A>G; p.D90G and c.H86G>T; p.V396F) in a com- 



pound heterozygous state in the affected siblings. The 
variants cosegregated with the disease in the entire family 
(Fig. 1) and were absent in 400 control chromosomes. 
Mutation D90G lay within the predicted equatorial 
domain and mutation V396F lay within the predicted 
intermediate domain which connects the equatorial and 



D90G 



V396F 



MKKS 



570 aa 



■ Equatorial domain □ Intermediate domain 
□ Apical domain □ Insertion 

g p.D90G 

1 

Human GGYVCTTSQSSALLSHLLVTHPILKILTASIQNHVSSFSDCGLFTAILCC 
Pan troglodytes GGYVCTTSQSSALLSHLLVTHPILKILTTSIQNHVSSFSDCGLFTAILCC 
Canis lupus f amiliaris GGCVCTTSQSSALLSNLPVTHPILKILTTSMQNHVSCFRDCGLFTAILCC 
Mus musculus GGCVYTTSQSSALLRNLSVTHPVLKILTSSVQNHVSCFSDCGLFTAILCC 
Gallus gallus GGYVCLTSQSSAILGHLSVSHPVLKVLTASVQNHVSRFSDCGLFTAILCC 
Xenopus tropicalis GGCVLTTSQSSALFNSFSVSKPVAKLLVASIRNHI SCFSDSGLFAASLCC 

Danio rerio GGHVLTTSTSTALLKRLEMSEPLLKLISTALQHHTTRYSDSGLFMGIFTL 

' 'p.'v396F' ' 

Human SLLLCNRNDTAWDELKLTCQTALHVLQLTLKEPWALLGGGCTETHLAAYI 
Pan troglodytes SLLLCNRNDTAWDELKLTCQTALHVLQLTLKEPWALLGGGCTETHLAAYI 
Canis lupus familiaris SLLLCNRNDTAWDELKLTCQTALHVLQLTIKEPCVLLGGGCTETHLAAYI 
Mus musculus TLLLCSRNDTAWEELKLTCQTAMHVLQLTIKEPWVLLGGGCTETHLAAYV 
Gallus gallus SLILCNRNETTWDELKRACETAEHVLQLTLKEPLALLGGGCTETHLASYL 
Xenopus tropicalis SFVLCNRNETSLKELMRTCEAAERVLQLTLKNPWVLLGGGCTETHLAALL 
Danio rerio TMVLCHRNETMLEELKMTCQRAEHVLRLTLREPYALLGGGCTETQLATHI 



Figure 3. MKKS protein depiction. (A) Schematic representation of the identified variants within the MKKS domains, including the typical 
chaperonin group II domains (equatorial, intermediate, and apical). The domain organization was modified from (Stoetzel et al. 2007). (B) 
Alignment of the orthologs from different species showing conservation of the mutated residues. An * (asterisk) indicates positions which have a 
single, fully conserved residue. A: (colon) indicates conservation between groups of strongly similar properties. A. (period) indicates conservation 
between groups of weakly similar properties. 
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the apical domains via flexible hinges (Stone et al. 2000) 
(Fig. 3A). Alignment of MKKS amino acid sequence of 
various orthologs showed that the substituted residues 
(aspartic acid at position 90 and valine at position 396) 
are highly conserved across species from different 
evolutionary branches (Fig. 3B). In silico tools predicted 
that both p.D90G and p.V396F are probably damaging by 
Polyphen-2 (score = 1 for p.D90G and score = 0.987 for 
p.V396F) and SIFT (score = 0.01 for p.D90G and 
score = 0 for p.V396F). 

Epistasis evaluation 

Whole-exome data were subjected to exhaustive evalua- 
tion paying special attention in the identification of addi- 
tional variants that may be acting as second-site modifiers 
making patient's 11:4 BBS phenotype more severe. 
Although several SNPs were identified in the remaining 
16 known BBS genes: BBS1 (MIM 209901), BBS2 (MIM 
606151), ARL6 (BBS3; MIM 608845), BBS4 (MIM 
600374), BBS5 (MIM 603650), MKKS (BBS6), BBS7 
(MIM 607590), TTC8 (BBS8; MIM 608132), BBS9 (MIM 
607968), BBS10 (MIM 610148), TRIM32 (BBS11; 
MIM 602290), BBS12 (MIM 610683), MKS1 (BBS 13; 
MIM 609883), CEP290 (BBS 14; MIM 610142), SDCCAG8 
(BBS 16; MIM 613524), LZTFL1 (BBS 17; MIM 606568), 
and INPP5E (MIM 613037) (Table SI); none of them was 
predicted to be pathogenic. 

To explore the possibility that the polycystic kidney 
disease (PKD), affecting only patient 11:4, may be 
caused by mutations in other loci such as PKD genes: 
PKD1 (MIM 601313), PKD2 (MIM 173910), and 
PKHD1 (MIM 606702); or nephronophthisis genes: 
NPHP1 (MIM 607100), INVS (NPHP2; MIM 243305), 
NPHP3 (MIM 608002), NPHP4 (MIM 607215), IQCB1 
(NPHP5; MIM 609237), GLIS2 (NPHP7; MIM 608539), 
RPGRIP1L (NPHP8; MIM 610937), NEK8 (NPHP9; 
MIM 609799), NEK1 (MIM 604588), MKS1 (MIM 
609883), and TMEM67 (NPHP11; MIM 609884) (Table 
SI), we evaluated the exome data and found one novel 
variant in intron 8-9 of the NPHP4 gene 
(c.992 + 71G>T). This variant has never been reported 
in databases such as dbSNP and EVS. As coverage 
information for this position was not clear in the EVS 
data set, we checked the absence of c.992 + 71G>T in 
additional genomic databases which also include intron- 
ic variants (1000 genomes and 5000 genomes). Family 
segregation showed the presence of the heterozygous 
NPHP4 variant only in patient 11:4 and in his unaf- 
fected parent (1:1) (Fig. 1). In silico tools predicted that 
the substitution of a G > T in this position activates a 
cryptic splice acceptor site in the intron 8-9 of NPHP4 
(Table 3). The complete coding sequence was scanned 



Table 3. Splice acceptor site prediction scores for NPHP4 
c.992 + 71G>T mutation versus wild type. 





Wild type 


c.992 + 71G>T 




aggtcatttgtg 


intronEXON 


Allele 


catgtcaggtgt 


aggtcatttgttcatgtcagGTGT 


NNSPI ICF 

1 \J 1 \ -J 1 1 1 v 1 


o 


0.71 


ESEfinder 






3SS_U2_human 


0 


7.735500 


(threshold: 6.632) 






3SS_U2_mouse 


0 


7.26420 


(threshold: 6.724) 






BranchSite 


0 


2.07090 (tgttcat) 


(threshold: 0) 






NetGene2 


0 


0.28 



Sequence variants are highlighted in bold; potential splice acceptor 
site are underlined; predicted exonic sequence is in capital letters. 
3SS_U2_Human: 3's splice sites (acceptor) of human (U2 type). 
3SS_U2_mouse: 3's splice sites (acceptor) of mouse (U2 type). Branch 
site: mammalian branch site (U2 type). NPHP4 Genbank accession 
number: NG_01 1724.2. 

to eliminate the possibility of the presence of other 
mutation in NPHP4. The absence in 400 control 
chromosomes and the results of the in silico prediction 
supported the pathogenic role of the NPHP4 c.992 + 
71G>T variant. 

Network analysis 

In order to explore the possible physical interactions 
among the BBS proteins and NPHP4, network enrich- 
ment analysis, using the program SNOW (Minguez et al. 
2009), has been performed. Figure 4 shows the proteins 
significantly connected (significantly smallest number of 
components, P- value = 0.01) allowing the introduction of 
one extra connecting node among the proteins studied. 

The network documents the dense network of physical 
protein-protein interactions that connect BBS proteins 
and also documents how NPHP4 connects to several BBS 
proteins through different intermediates. Of special inter- 
est is RPGR, which connects NPHP4 to CEP290. RPGR, 
the RP GTPase regulator, is also involved in RP X-Linked, 
cone-rod dystrophy X-linked, and macular dystrophy 
X-linked. Another interesting gene is INVS, Inversin, 
which encodes a protein that protein may function in 
renal tubular development and function. Such connec- 
tions reinforce the possible role of NPHP4 as modulator 
of the penetrance of the disease. 

Discussion 

In this report, a Spanish BBS family with three affected sib- 
lings is described. The mode of inheritance and the main 
clinical features correspond to autosomal recessive BBS. 
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Exome analysis of four individuals (1:1, 1:2, 11:3, and 
11:4) led to the identification of two novel compound het- 
erozygous mutations in the MKKS gene (c.269A>G; 
p.D90G and c.H86G>T; p.V396F) of affected individuals. 
These variants were absent in 200 control individuals and 
showed segregation with disease in the entire family. 
Mutation D90G lay within the predicted equatorial 
domain, the most conserved region among group I and II 
chaperonins. In the group II chaperonin family, this 
domain is responsible for ATP hydrolysis and the substi- 
tution might alter its structure. Mutation V396F lay 
within the predicted intermediate domain which connects 
the equatorial and the apical domains via flexible hinges 
(Stone et al. 2000). Bioinformatic analysis predicted path- 
ogenic consequences for both missense mutations. 

The molecular diagnosis is generally helpful to confirm 
a clinical diagnosis. Although many studies have dissected 
the clinical overlaps due to MKKS mutations, genotype- 
phenotype correlations are still not well understood. The 
clinical evaluation of the affected members of the family 
did not reveal hydrometrocolpos, suggesting that these 
two missense mutations in the MKKS gene did not cause 
genital malformations typical of MKKS phenotype. In this 
family, mutations in MKKS resulted in a spectrum of BBS 
phenotypes, ranging from BBS with severe renal involve- 
ment in patient 11:4 to milder forms of BBS (patients 11:2 
and 11:5). The expression of the phenotype and the dis- 
ease progression can vary greatly from patient to patient, 
even among members of the same family. Whereas retinal 
dystrophy, digit anomaly, and learning disabilities were 



highly penetrant traits for all patients, kidney abnormali- 
ties and obesity showed incomplete penetrance in our 
family. 

The index patient investigated here had previously 
undergone selected genotyping (APEX analysis, Asper 
http://www.asperbio.com/asper-ophthalmics), but this 
mutational screening approach failed to identify the 
underlying gene defect in this family. The implementation 
of a reliable diagnostic system able to detect novel dis- 
ease-causing mutations, even in genes not previously 
associated with an assumed clinical diagnosis, is necessary. 
Exome sequencing has proven to be an important diag- 
nostic tool for disorders that are characterized by signifi- 
cant genetic heterogeneity. In the case of BBS, besides the 
high number of genes involved, oligogenic inheritance is 
also well documented (Katsanis 2004) adding a layer of 
complexity to the genetic characterization of such 
patients. 

One advantage of analyzing the whole exome simulta- 
neously is the unique possibility to investigate the 
involvement of second-site mutations which may be mod- 
ulating the expression of the BBS phenotype (Katsanis 
et al. 2001; Badano et al. 2003; Beales et al. 2003; Katsanis 
2004; Hjortshoj et al. 2010). However, the use of unbi- 
ased methodologies often produces many candidates that 
must be filtered out with in silico prioritization tech- 
niques. Thus, network analysis captures the relationships 
of mutated genes with already known disease genes allow- 
ing a rational prioritization of candidate genes. The com- 
bination of NGS with network analysis approach 
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conducted here allowed us to identify a novel variant - 
located on intron 8-9 of NPHP4 - which was predicted 
to enhance the use of a cryptic splice acceptor site. This 
would probably cause the introduction of a premature 
termination codon and the reduction in NPHP4 mRNA 
levels. Only the most severely affected patient (11:4) car- 
ries the variant (c.992 + 71G>T) in heterozygous state. In 
addition, the absence of the variant in control individuals 
seems to indicate that it is a pathogenic allele. 

NPHP4 is a ciliary protein that is known to belong to 
a multifunctional complex and colocalizes with RPGPJP1, 
RPGR, and the serologically defined colon cancer antigen- 
8 (SDCCAG8), a protein thought to partake in the 
RPGRIP1 interactome and implicated also in retinal-renal 
ciliopathies (Roepman et al. 2005; Schaefer et al. 2011; 
Won et al. 2011; Patil et al. 2012). Mutations in NPHP4 
have been associated with related ciliopathies such as 
NPH (Mollet et al. 2002; Hoefele et al. 2004) and SLS 
(Otto et al. 2002; Schuermann et al. 2002). 

Few splice site variants have been found to modulate 
the pathogenesis of BBS. One example is the allele C430T 
of MGC1203 which enhances the use of a cryptic splice 
acceptor site, causing the introduction of a premature ter- 
mination codon and the reduction in the steady-state 
MGC1203 mRNA levels (Badano et al. 2006b). Although 
the MGC1203 mutations are probably insufficient to 
cause BBS, this gene may be involved in the pathogenesis 
of BBS by contributing hypomorphic mutations to an 
already sensitized genetic background. The authors con- 
clude that suppression of MGC1203 exerts an epistatic 
effect on the developmental BBS phenotype. It is tempt- 
ing to speculate that NPHP4 allele reported herein may 
also show a similar effect. We hypothesize that the 
NPHP4 variant may be acting as a second-site modifier 
affecting the likelihood of renal degeneration in the con- 
text of the other MKKS mutations. This hypothesis is also 
supported by the known physical relationships revealed 
by network analysis. The MKKS protein forms a signifi- 
cant cluster of proteins with several BBS proteins, as 
shown in Figure 4. The severe BBS phenotype observed in 
patient 11:4 may be due to an epistatic effect of MKKS 
and NPHP4 mutations. This would provide an explana- 
tion for one of the distinct clinical manifestations 
observed in this family. Moreover, it is possible that other 
intrafamilial phenotypic variations may be due to altera- 
tions in additional loci as well as nonstrictly genetic fac- 
tors not evaluated in this study. Therefore, our findings 
warrant further studies to evaluate the role of putative 
disease-modifying variants. 

Identifying the presence of epistatic variants can be chal- 
lenging, given that a large number of common and rare 
alleles - located in known and in novel genes - has been 
shown to influence the expressivity of the BBS phenotype. 



Our data demonstrate that NGS is a powerful approach to 
identify rare and high penetrant disease variants and to 
study the genetic contribution to phenotypic variability. 

This genetic diagnostic tool can now be applied to large 
cohorts of Mendelian and oligogenic disorders and should 
rapidly provide the molecular diagnosis and the preva- 
lence of associated genes. 

In summary, we show how WES has clearly improved 
the molecular diagnosis of heterogeneous disorders such 
as BBS. In addition, our study highlights the usefulness of 
NGS approaches for the dissection of epistatic phenom- 
ena and provides promising findings to decipher the 
genetic basis of phenotypic variability. 
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